#include "Resolver.h"
#include <cmath>

Resolver::Resolver(ifstream &datos)
{
	// Orden dentro del archivo de entrada de los valores:
	// Radio B
	// Radio A
	// DeltaTeta
	// DeltaRadio
	// Vector
	std::vector<double> tempBorde;
	datos >> _radioB;
	datos >> _radioA;
	datos >> _delta_theta;
	datos >> _delta_radio;
	int puntosExt = 360 / (int)_delta_theta;
	double temp;
	for (int i = 0; i < puntosExt; i++)
	{
		datos >> temp;
		tempBorde.push_back(temp);
	}
	int cantPuntos = (int)(_radioA / _delta_radio);
	cantPuntos *= puntosExt;

	// Ahora tengo en cantPuntos la cantidad total de los puntos de la matriz.
	double **matriz = new double*[cantPuntos];
	for (int i = 0; i < cantPuntos; i++)
		matriz[i] = new double[cantPuntos+1];
	// Utilizamos cantPuntos + 1 porque asi nos queda una colmuna para poder guardar los resultados. FG
	// LLenamos de ceros la matriz
	for(int i = 0; i < cantPuntos; i++)
	{
		for(int j = 0; j < cantPuntos + 1; j++)
		{
			matriz[i][j] = 0;
		}
	}
	// Armamos la matriz
	int radio = 0;
	int particion_angulos = 360 / (int)_delta_theta;
	for(int fila = 0; fila < cantPuntos; fila++)
	{
		radio = (int)((fila / particion_angulos) * _delta_radio + _radioB);
		matriz[fila][fila] = 1.0;
		if (fila < particion_angulos)
		    matriz[fila][cantPuntos] = 1500.0;

		else if (fila + particion_angulos > cantPuntos - 1)
				matriz[fila][cantPuntos] = tempBorde[fila % particion_angulos];

		else
		{
		    //caso angulo == 0
		    if(fila % particion_angulos == 0)
		    {
			// Como es el primero, el anterior es el ultimo
		        matriz[fila][fila + particion_angulos - 1] = valor_k_masmenos_uno(radio);
		        matriz[fila][fila+1] = valor_k_masmenos_uno(radio);
		        if( (fila - particion_angulos) < particion_angulos )
		        {
		            // Me encuentro en el primer radio
		            matriz[fila][cantPuntos] -= 1500.0 * valor_j_menos_uno(radio); // tengo un dato que paso para el otro lado
		            if (fila + particion_angulos >= (cantPuntos - particion_angulos) )
		             // Si tengo el caso de radio máximo
		                matriz[fila][cantPuntos] -= tempBorde[fila % particion_angulos] * valor_j_mas_uno(radio); // también
		            else
		                matriz[fila][fila + particion_angulos] = valor_j_mas_uno(radio);
		            
		        }
		        else if (fila + particion_angulos >= (cantPuntos - particion_angulos)) 
                // si tengo el caso de radio máximo
		        {
		            matriz[fila][cantPuntos] -= tempBorde[fila % particion_angulos] * valor_j_mas_uno(radio); // también
		            matriz[fila][fila - particion_angulos] = valor_j_menos_uno(radio);
		        }
		        else
		        {
		            matriz[fila][fila - particion_angulos] = valor_j_menos_uno(radio);
		            matriz[fila][fila + particion_angulos] = valor_j_mas_uno(radio);
		        }
		        // El gran bloque de if está porque estamos analizando varios casos disjuntos. 
			// Si no se da ninguna de las posibilidades, entonces ponemos los valores 
			// conocidos
		    }
		    //caso angulo máximo
		    else if((fila % particion_angulos) == particion_angulos - 1)
		    {
			// Tamos viendo si es el ultimo angulo del radio.
		        matriz[fila][fila + 1 - particion_angulos] = valor_k_masmenos_uno(radio);
		        matriz[fila][fila-1] = valor_k_masmenos_uno(radio);
		        if( (fila - particion_angulos) < particion_angulos) // si me encuentro con el radio 0
		        {
		            matriz[fila][cantPuntos] -= 1500.0 * valor_j_menos_uno(radio); // tengo un dato que paso para el otro lado
		            if (fila + particion_angulos >= (cantPuntos - particion_angulos)) // si tengo el caso de radio máximo
		                matriz[fila][cantPuntos] -= tempBorde[fila % particion_angulos] * valor_j_mas_uno(radio); // también
		            else
		                matriz[fila][fila + particion_angulos] = valor_j_mas_uno(radio);
		        }
		        else if (fila + particion_angulos >= (cantPuntos - particion_angulos)) // si tengo el caso de radio máximo
		        {
		            matriz[fila][cantPuntos] -= tempBorde[fila % particion_angulos] * valor_j_mas_uno(radio); // también
		            matriz[fila][fila - particion_angulos] = valor_j_menos_uno(radio);
		        }
		        else
		        {
		            matriz[fila][fila - particion_angulos] = valor_j_menos_uno(radio);
		            matriz[fila][fila + particion_angulos] = valor_j_mas_uno(radio);
		        }
		    }
		    else
		    {
		        if( (fila - particion_angulos) < particion_angulos) // si me encuentro con el radio 1
		        {
		            matriz[fila][cantPuntos] -= 1500.0  * valor_j_menos_uno(radio); // tengo un dato que paso para el otro lado
		            if (fila + particion_angulos >= (cantPuntos - particion_angulos)) // si tengo el caso de radio máximo
		                matriz[fila][cantPuntos] -= tempBorde[fila % particion_angulos] * valor_j_mas_uno(radio); // también
		            else
		                matriz[fila][fila + particion_angulos] = valor_j_mas_uno(radio);
		        }
		        else if (fila + particion_angulos >= (cantPuntos - particion_angulos)) // si tengo el caso de radio máximo
		        {
		            matriz[fila][cantPuntos] -= tempBorde[fila % particion_angulos] * valor_j_mas_uno(radio); // también
		            matriz[fila][fila - particion_angulos] = valor_j_menos_uno(radio);
		        }
		        else
		        {
		            matriz[fila][fila - particion_angulos] = valor_j_menos_uno(radio);
		            matriz[fila][fila + particion_angulos] = valor_j_mas_uno(radio);
		        }

		        // Caso normal de un punto en el medio
		        matriz[fila][fila-1] = valor_k_masmenos_uno(radio);
		        matriz[fila][fila+1] = valor_k_masmenos_uno(radio);
		    }
		} 
	}
	Matriz ret(matriz, cantPuntos, cantPuntos+1);
	_matrix = ret;
}

void Resolver::terminarTp()
{
    ofstream matrizss("matriz.txt");
    matrizss << _matrix;
    matrizss.close();
    
	_matrix.triangular();
	
	ofstream matrizss_t("matriz_triang.txt");
    matrizss_t << _matrix;
    matrizss_t.close();
    
	_resultados = _matrix.resolver();
	
	ofstream matrizss_r("matriz_resol.txt");
    matrizss_r << _matrix;
    matrizss_r.close();

	Matriz matFinal = generarMatrizDiscretizada();
	escritura(matFinal);
}

void Resolver::escritura(Matriz& matriz)
{
	ofstream file("salida.txt");

	std::vector<Coordenada> isoterma = _matrix.calcularIsoterma(_delta_theta, _delta_radio, _radioB);
	Matriz mat = generarMatrizDiscretizada();
	
	file << "Solución del trabajo práctico número 2" << endl;
	file << "Para los parámetros" << endl;
	file << "Delta theta = " << _delta_theta << endl;
	file << "Delta r = " << _delta_radio << endl;
	file << "Radio de la pared interna del horno = " << _radioB <<endl;
	file << "La matriz de temperatura es: " << endl << endl;
	file << mat << endl;
	file << "La función isoterma pasa por los puntos (cartesianos)" << endl;
	double x, y;
	for(int k = 0; k < isoterma.size(); k++)
	{
		x = (isoterma[k].radio + _radioB)*cos(isoterma[k].angulo * (PI / 180.0));
		y = (isoterma[k].radio + _radioB)*sin(isoterma[k].angulo * (PI / 180.0));
		file << "(" << x << " , " << y << ")" << endl;
	}
	file.close();

	imprimir_matlab();
}
	


double Resolver::valor_k_masmenos_uno(double r) 
{
	double numerador, denominador;
	numerador = pow(_delta_radio, 2);
	denominador = r*_delta_radio - 2.0*pow(r, 2);
	denominador *= pow(_delta_theta, 2);
	denominador -= 2.0*_delta_radio;
	return numerador / denominador;
}

double Resolver::valor_j_mas_uno(double r)
{
    double numerador, denominador;
	numerador = pow(r, 2) * pow(_delta_theta, 2);
	denominador = r*_delta_radio - 2.0*pow(r, 2);
	denominador *= pow(_delta_theta, 2);
	denominador -= 2.0*_delta_radio;
	return numerador / denominador;
}

double Resolver::valor_j_menos_uno(double r)
{
	double numerador, denominador;
	numerador = pow(r, 2) * pow(_delta_theta, 2);
	numerador -= r*_delta_radio*pow(_delta_theta, 2);
	denominador = r*_delta_radio - 2.0*pow(r, 2);
	denominador *= pow(_delta_theta, 2);
	denominador -= 2.0*_delta_radio;
	return numerador / denominador;
}

Matriz Resolver::generarMatrizDiscretizada()
{
	int cant_radios = (int) (_radioA / _delta_radio);
	int cant_grados = (int) (360.0 / _delta_theta);
	double **matriz_discreta = new double*[cant_radios];
	for(int i = 0; i < cant_radios; i++)
		matriz_discreta[i] = new double[cant_grados];

	for(int radios = 0; radios < cant_radios; radios++)
		for(int grados = 0; grados < cant_grados; grados++)
			matriz_discreta[radios][grados] = _resultados[radios*cant_grados+grados];

	return Matriz(matriz_discreta, cant_radios, cant_grados);
}

void Resolver::imprimir_matlab()
{
	ofstream matfile("tp2.m");
	Matriz mat = generarMatrizDiscretizada();
	imprimir_matriz_x(matfile, mat);
	imprimir_matriz_y(matfile, mat);
	imprimir_matriz_z(matfile, mat);
	matfile << "surf(X, Y, Z);" << endl;
	matfile.close();
}

void Resolver::imprimir_matriz_x(ofstream &of, Matriz &mat)
{
	of << "X = [";
	double angulo, radio, x;
	for(int i = 0; i < mat.cantFilas(); i++)
	{
		radio = i*_delta_radio + _radioB;
		for(int j = 0; j < mat.cantColumnas(); j++)
		{
			angulo = j*_delta_theta;
			angulo *= (PI / 180.0);
			x = radio*cos(angulo);
			of << x << "\t";
		}
		// se agrega nuevamente el primer punto de la discretización por la circularidad del área
		angulo = 2.0*PI - 0.0000001;
		x = radio*cos(angulo);
        of << x << "\t";
		of << ";\t";
	}
	of << "];" << endl;
}

void Resolver::imprimir_matriz_y(ofstream &of, Matriz &mat)
{
	of << "Y = [";
	double angulo, radio, y;
	for(int i = 0; i < mat.cantFilas(); i++)
	{
		radio = i*_delta_radio + _radioB;
		for(int j = 0; j < mat.cantColumnas(); j++)
		{
			angulo = j*_delta_theta;
			angulo *= (PI / 180.0);
			y = radio*sin(angulo);
			of << y << "\t";
		}
		// se agrega nuevamente el primer punto de la discretización por la circularidad del área
		angulo = 2.0*PI - 0.0000001;
		y = radio*sin(angulo);
		of << y << "\t";
		of << ";\t";
	}
	of << "];" << endl;
}

void Resolver::imprimir_matriz_z(ofstream &of, Matriz &mat)
{
	of << "Z = [";
	Matriz m = generarMatrizDiscretizada();
	for(int i = 0; i < m.cantFilas(); i++)
	{
		for(int j = 0; j < m.cantColumnas(); j++)
			of << m.valor(i,j) << "\t";
		of << m.valor(i, 0) << ";\t";
	}
	of << "];" << endl;
}
